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, The purpose of this paper is twofold: from one side we provide a 

Pi I general survey to the viscoelastic models constructed via fractional 

calculus and from the other side we intend to analyze the basic 
fractional models as far as their creep, relaxation and viscosity 
properties are considered. The basic models are those that generalize 
' via derivatives of fractional order the classical mechanical models 

. characterized by two, three and four parameters, that we refer to 

j. I as Kelvin-Voigt, Maxwell, Zener, anti-Zener and Burgers. For each 

' fractional model we provide plots of the creep compliance, relaxation 

modulus and effective viscosity in non dimensional form in terms of 
a suitable time scale for different values of the order of fractional 
derivative. We also discuss the role of the order of fractional derivative 
in modifying the properties of the classical models. 
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1 Introduction 



A topic of continuum mechanics, where fractional calculus is suited to be 
applied, is without doubt the linear theory of viscoelasticity. In fact, an 
increasing number of authors have used fractional calculus as an empirical 
method of describing the properties of linear viscoelastic materials. A wide 
bibliography up to nowadays is contained in the recent book by Mainardi 
[27j including an historical perspective up to 1980's. 

The purpose of this paper is to provide, after a general survey to the linear 
theory of viscoelasticity, a (more) systematic discussion and a graphical 
representation of the main properties of the basic models described by stress- 
strain relationships of fractional order. The properties under discussion 
concern the standard creep and relaxation tests that have a relevance in 
experiments. 

The plan of paper is as following. In Section 2 we recall the essential notions 
of linear viscoelasticity in order to present our notations for the analog 
mechanical models. We limit our attention to the basic mechanical models, 
characterized by two, three and four parameters, that we refer to as Kelvin- 
Voigt, Maxwell, Zener, anti-Zener and Burgers. 

In Section 3 we consider our main topic concerning the creep, relaxation and 
viscosity properties of the previous basic models generalized by replacing in 
their differential constitutive equations the derivatives of integer order 1 and 
2 with derivatives of fractional order v and 1 + z/ respectively, with < z/ < 1. 
We provide the analytical expressions and the plots of the creep compliance, 
relaxation modulus and effective viscosity for all the considered fractional 
models. 

We also enclose two Appendices for providing the readers with the essential 
notions of fractional derivative and with a discussion on initial conditions. 

2 Essentials of linear viscoelasticity 

In this section we present the fundamentals of linear viscoelasticity restricting 
our attention to the one-axial case and assuming that the viscoelastic body 
is quiescent for all times prior to some starting instant that we assume as 
t = 0. 
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For the sake of convenience both stress a{t) and strain e{t) are intended to 
be normahzed, i.e. scaled with respect to a suitable reference state {do , eo} . 

2.1 Generalities 

According to the linear theory, the viscoelastic body can be considered as a 
linear system with the stress (or strain) as the excitation function (input) and 
the strain (or stress) as the response function (output). In this respect, the 
response functions to an excitation expressed by the Heaviside step function 
Q{t) are known to play a fundamental role both from a mathematical and 
physical point of view. We denote by J{t) the strain response to the unit 
step of stress, according to the creep test and by G(t) the stress response to 
a unit step of strain, according to the relaxation test. 

The functions J{t) and G{t) are usually referred to as the creep compliance 
and relaxation modulus respectively, or, simply, the material functions of the 
viscoelastic body. In view of the causality requirement, both functions are 
vanishing for t < 0. 

The limiting values of the material functions for t — > 0^ and t — )■ +oo are 
related to the instantaneous (or glass) and equilibrium behaviours of the 
viscoelastic body, respectively. As a consequence, it is usual to denote Jg : = 
J(0+) the glass compliance, Jg := J(+oo) the equilibrium compliance, and 
Gg := G'(O^) the glass modulus Ge '.= G(+oo) the equilibrium modulus. As a 
matter of fact, both the material functions are non-negative. Furthermore, 
for < t < +00 , J{t) is a non decreasing function and G{t) is a non 
increasing function. 

The monotonicity properties of J{t) and G{t) are related respectively to 
the physical phenomena of strain creep and stress relaxation. We also note 
that in some cases the material functions can contain terms represented by 
generalized functions (distributions) in the sense of Gel'fand-Shilov [12] or 
pseudo-functions in the sense of Doetsch |10j . 

Under the hypotheses of causal histories, we get the stress-strain relation- 
ships 




(2.1) 
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where the passage to the RHS is justified if different iabihty is assumed 
for the stress-strain histories, see also the excellent book by Pipkin [33]. 
Being of convolution type, equations (2.1) can be conveniently treated by 
the technique of Laplace transforms so they read in the Laplace domain 



e(s) = s J{s) a{s) , cr(s) = s G{s) e{s) 
from which we derive the reciprocity relation 

1 

sJ(s) = ^ — . 
sG{s) 



(2.2) 



(2.3) 



Because of the limiting theorems for the Laplace transform, we deduce 
that Jg = 1/Gg, Je = 1/Ge, with the convention that and +oo are 
reciprocal to each other. The above remarkable relations allow us to classify 
the viscoelastic bodies according to their instantaneous and equilibrium 
responses in four types as stated by Caputo and Mainardi in their 1971 
review paper [S], see Table 2.1. 



Type 




Je 


Gg 


Ge 


I 


> 


< 00 


< 00 


> 


II 


> 


= 00 


< 00 


= 


III 


= 


< 00 


= 00 


> 


IV 


= 


= 00 


= 00 


= 



Table 2.1 The four types of viscoelasticity. 



We note that the viscoelastic bodies of type I exhibit both instantaneous 
and equilibrium elasticity, so their behaviour appears close to the purely 
elastic one for sufficiently short and long times. The bodies of type II and 
IV exhibit a complete stress relaxation (at constant strain) since Ge = and 
an infinite strain creep (at constant stress) since Jg = cxd , so they do not 
present equilibrium elasticity. Finally, the bodies of type III and IV do not 
present instantaneous elasticity since Jg = {Gg = 00). Other properties 
will be pointed out later on. 



2.2 The mechanical models 

To get some feeling for linear viscoelastic behaviour, it is useful to consider 
the simpler behaviour of analog mechanical models. They are constructed 
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from linear springs and dashpots, disposed singly and in branches of two 
(in series or in parallel). As analog of stress and strain, we use the total 
extending force and the total extension. We note that when two elements 
are combined in series [in parallel], their compliances [moduli] are additive. 
This can be stated as a combination rule: creep compliances add in series, 
while relaxation moduli add in parallel. 

The mechanical models play an important role in the literature which is 
justified by the historical development. In fact, the early theories were 
established with the aid of these models, which are still helpful to visualize 
properties and laws of the general theory, using the combination rule. 

Now, it is worthwhile to consider the simple models of Figs. 1, 2, 3 
providing their governing stress-strain relations along with the related 
material functions. 



The Hooke and the Newton models. The spring (Fig. la) is the elastic 
(or storage) element, as for it the force is proportional to the extension; it 
represents a perfect elastic body obeying the Hooke law (ideal sohd). This 
model is thus referred to as the Hooke model. If we denote by m the pertinent 
elastic modulus we have 

Hooke model : a{t)=me{t), and I 'i^jl 2^^^' (2-4) 



G{t) =m. 

In this case we have no creep and no relaxation so the creep compliance 
and the relaxation modulus are constant functions: J{t) = Jg = — 1/m; 
G{t) = Gg = Ge = m. 

The dashpot (Fig. lb) is the viscous (or dissipative) element, the force 
being proportional to rate of extension; it represents a perfectly viscous body 
obeying the Newton law (perfect liquid). This model is thus referred to as 
the Newton model. If we denote by bi the pertinent viscosity coefficient, we 
have 

..^ / Jit) = t/bi , 



Newton model : cr(t) = 6i — , and < / (2.5) 

^ ^ dt \ G{t) = bi5{t) . ^ ^ 

In this case we have a linear creep J{t) — J+t and instantaneous relaxation 
G{t) = (?_ S{t) with G_ = 1/J+ = 6i. 

We note that the Hooke and Newton models represent the limiting cases of 
viscoelastic bodies of type / and IV, respectively. 
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The Kelvin— Voigt and the Maxwell models. A branch constituted by 

a spring in parallel with a dashpot is known as the Kelvin-Voigt model (Fig. 
Ic). We have 



Kelvin — Voigt model : a{t) = m €{t) + bi ^ , 

CLL 



(2.6a) 



and 



J{t) = Ji 



1 - e-V^. 



T ^ ^1 

, Jl = — , — 

m m 



(2.66) 



G{t) ^Ge + G^ S{t) , Ge^m, G-^bi, 
where is referred to as the retardation time. 

A branch constituted by a spring in series with a dashpot is known as the 
Maa;w;e// model (Fig. Id). We have 



, , , , da de 
Maxwell model : alt) + ai — = Oi — , 
^ ^ dt dt' 



(2.7a) 



and 



, , ai 1 

■^(^J — Jg + J+t , "^9 ^ TT ' + ^ 77 

Oi Oi 



G{t) = Gie' 



ai 



where is referred to as the relaxation time. 



(2.76) 



The Voigt and the Maxwell models are thus the simplest viscoelastic bodies 
of type /// and //, respectively. The Voigt model exhibits an exponential 
(reversible) strain creep but no stress relaxation; it is also referred to as the 
retardation element. The Maxwell model exhibits an exponential (reversible) 
stress relaxation and a linear (non reversible) strain creep; it is also referred 
to as the relaxation element. 

The Zener and the anti— Zener models. Based on the combination rule, 
we can continue the previous procedure in order to construct the simplest 
models of type / and IV that require three parameters. 

The simplest viscoelastic body of type / is obtained by adding a spring either 
in series to a Voigt model or in parallel to a Maxwell model (Fig. 2a and 
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Figure 1: The representations of the basic mechanical models: a) spring for 
Hooke, b) dashpot for Newton, c) spring and dashpot in parallel for Voigt, 
d) spring and dashpot in series for Maxwell. 



Fig. 2b, respectively). So doing, according to the combination rule, we 
add a positive constant both to the Voigt-like creep compliance and to the 
Maxwell-like relaxation modulus so that we obtain Jg > and Ge > ■ Such 



a model was introduced by Zener 
Linear Solid {S.L.S.). We have 



2l with the denomination of Standard 



Zener model 



l + ai 



d 
di 



a{t) 



m + hi 



d 
di 



e{t) 



(2.8a) 



and 

G{t) = Ge + Gie"^/^-, 



1 - e-^/^^ 



Jo 



Ji 



m, Gi 



m 

h 

ai 



hi 



m' 
ai 



(2.86) 



We point out the condition < m < 6i/ai in order Ji,Gi be positive and 
hence < < Jg < oo and < Gg < < oo . As a consequence, we note 
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that, for the S.L.S. model, the retardation time must be greater than the 
relaxation time, i.e. < < < oo . 




Figure 2: The mechanical representations of the Zener model, see a), b) and 
of the anti-Zener model, see c), d), where: a) spring in series with Voigt, b) 
spring in parallel with Maxwell, c) dashpot in series with Voigt, d) dashpot 
in parallel with Maxwell. 



Also the simplest viscoelastic body of type IV requires three parameters, 
i. e. ai , &i , &2 ; it is obtained adding a dashpot either in series to a Voigt 
model or in parallel to a Maxwell model (Fig. 2c and Fig 2d, respectively). 
According to the combination rule, we add a linear term to the Voigt-like 
creep compliance and a delta impulsive term to the Maxwell-like relaxation 
modulus so that we obtain Jg = oo and Gg — oo . We may refer to this model 
to as the anti-Zener model. We have 



anti — Zener model 
and 

J(i) = J+i + Ji 

G{t)=G_5{t) + Gie 



1 + ai 



d 
~dt 



a{t) 



e{t) , (2.9a) 



1 J _0'l h 

bi hi b{ 

U_ — — , — 2) ^o- 



b_2 

bl' 



(2.96) 



ai. 
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We point out the condition < < in order Ji,Gi be positive. As 

a consequence, we note that, for the anti-Zener model, the relaxation time 
must be greater than the retardation time, i.e. < < < 00 , on the 
contrary of the Zener [S.L.S.) model. 

In Fig. 2 we exhibit the mechanical representations of the Zener model 
(2.8a)-(2.8b), see a), b), and of the anti-Zener model (2.9a)-(2.9b), see c), 
d). 

Remark: We note that the constitutive equation of the anti-Zener model 
is formally obtained from that of the Zener model by replacing the strain 
e{t) by the stain-rate e{t). However the Zener model, introduced by Zener 
in 1948 |l2] for anelastic metals, was formerly introduced by Jeffreys with 
respect to bodily imperfection of elasticity in tidal friction, since from the 
first 1924 edition of his treatise on the Earth [T71 dB] , Sir Harold Jeffreys was 
usual to refer to the rheology of the Kelvin- Voigt model (that was suggested 
to him by Sir J. Larmor) to as firmoviscosity and to the rheology of the 
Maxwell model to as elastoviscosity. We observe that in the literature of 
rheology of viscoelastic fluids (including polymeric liquids) our anti-Zener 
model is (surprisingly for us) known as Jeffreys fluid, see e.g. the review 
paper by Bird and Wiest [2]. Presumably this is due to the replacement 
of the strain with the strain-rate (suitable for fluids) in the stress-strain 
relationship introduced by Jeffreys. As a matter of fact, in Earth rheology 
the Jeffreys model is known to be the creep model introduced by him in 
1958, see [IS], as generalization of the Lomnitz logarithmic creep law and 
well described in the subsequent editions of Jeffreys' treatise on the Earth. 
In view of above considerations, we are tempted to call our anti-Zener model 
as Standard Linear Fluid in analogy with the terminology Standar Linear 
Solid commonly adopted for the Zener model. 

The Burgers model. In Rheology literature it is customary to consider 
the so-called Burgers model, which is obtained by adding a dashpot or 
a spring to the representations of the Zener or of the anti-Zener model, 
respectively. Assuming the creep representation the dashpot or the spring 
is added in series, so the Burgers model results in a series combination 
of a Maxwell element with a Voigt element. Assuming the relaxation 
representation, the dashpot or the spring is added in parallel, so the Burgers 
model results in two Maxwell elements disposed in parallel. We refer the 
reader to Fig. 3 for the two mechanical representations of the Burgers model. 
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According to our general classification, the Burgers model is thus a four- 
element model of type II, defined by the four parameters {ai, 02, 61, 62}- 
We have 



Burgers model 



a{t)- 



d 



d^ 



e{t), (2.10a) 



so 



J{t) ^ Jg + J+t + Ji(l- e-^/^^) , 
G{t) = Gi e-*/^<^.i + G2 e-^/^^.2 . 



(2.106) 



We leave to the reader to express as an exercise the physical quantities Jg, 
J+, and Gi, r^.i, G2, Ta^2, in terms of the four parameters {ai, a2, 61, 62} in 
the operator equation (2.10a). 



The operator equation for the mechanical models Based on the 
combination rule, we can construct models whose material functions are of 
the following type 



J{t) =J, + 5^J„[l-e-^/^^.-]+J+i, 

n 

G{t) =Ge + ^G„e-^/^-.»+G_5(i), 



(2.11) 



where all the coefficient are non-negative, and interrelated because of the 
reciprocity relation (2.3) in the Laplace domain. We note that the four types 
of viscoelasticity of Table 2.1 are obtained from Eqs. (2.11) by taking into 
account that 



Je < 00 
Gg < 00 



J. = 



Jp = 00 



G_ = 0, Gg = 00 



J+^0, 
G_^Q. 



Appealing to the theory of Laplace transforms, we write 



+ 



G. 



sG{s) =(Ge + ^)- V-^^ + G_s, 

^ 1 + S T„_r,. 



(2.12) 



(2.13) 
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Figure 3: The mechanical representations of the Burgers modeh the creep 
representation (top) and the relaxation representation (bottom). 

where we have put (3 = J2n ^« • 

Furthermore, as a consequence of (2.13), sJ{s) and sG{s) turn out to be 
rational functions in C with simple poles and zeros on the negative real axis 
and, possibly, with a simple pole or with a simple zero at s = , respectively. 

In these cases the integral constitutive equations (2.1) can be written in 
differential form. Following Bland [3] with our notations, we obtain for these 
models 




(2.14) 



where q and p are integers with q = p or q = p+ l and m,ak,bk are non- 
negative constants, subjected to proper restrictions in order to meet the 



physical requirements of realizability. The general Eq. (2.14) is referred to 
as the operator equation of the mechanical models. 

In the Laplace domain, we thus get 



sJ{s) 



1 Pis] 



Pis) = 1 + J2(^ks' 



k=l 



sGis) Qis) 



where < 



(2.15) 



Qis) = m + ^6fcs'=. 



k=l 



with m > and q = p or q = p+ l. The polynomials at the numerator 
and denominator turn out to be Hurwitz polynomials (since they have no 
zeros for Re {s} > 0) whose zeros are alternating on the negative real axis 
is < 0). The least zero in absolute magnitude is a zero of Qis). The four 
types of viscoelasticity then correspond to whether the least zero is ( J+ 7^ 0) 
or is not ( J+ = 0) equal to zero and to whether the greatest zero in absolute 
magnitude is a zero of -P(s) (J^ 7^ 0) or a zero of Qis) (J^ = 0). 

In Table 2.2 we summarize the four cases, which are expected to occur in the 
operator equation (2.14), corresponding to the four types of viscoelasticity. 



Type 


Order 


m 


J9 






G. 


I 


q = p 


> 




m 








II 


q = p 


= 


ap/hp 





1/61 





III 


q = p + 1 


> 





m 





bq/ap 


IV 


q = p + 1 


= 








1/61 


hq/ap 



Table 2.2: The four cases of the operator equation. 

We recognize that for p = 1, Eq. (2.14) includes the operator equations for 
the classical models with two parameters: Voigt and Maxwell, illustrated in 
Fig. 1, and with three parameters: Zener and anti-Zener, illustrated in Fig. 
2. In fact we recover the Voigt model (type III) for m > and p = 0, g = 1, 
the Maxwell model (type II) for m = and p = g = 1, the Zener model 
(type I) for m > and p = g = 1, and the anti-Zener model (type IV) for 
m = and p = l,g = 2. Finally, with four parameters we can construct 
two models, the former with m = and p = g = 2, the latter with m > 
and p = 1, g = 2, referred by Bland [3] to as four-element models of the first 
kind and of the second kind, respectively. According to our convention they 
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are of type II and III, respectively. We have restricted our attention to the 
former model, the Burgers model (type II), illustrated in Fig. 3, because it 
has found numerous applications, specially in geo-sciences, see e.g. the books 
by Klausner |21j and by Carcione [9]. 



2.3 Complex modulus, effective modulus and effective 
viscosity 

In Earth rheology and seismology it is customary to write the one dimensional 
stress-strain relation in the Laplace domain in terms of a complex shear 
modulus /i(s) as 

a{s)=2il{s)eis), (2.16) 

that is expected to generalize the relation for a perfect elastic solid {Hooke 
model) 

a(t) = 2/ioe(t), (2.17) 
where /zq denotes the shear modulus. 

Adopting this notation we note comparing (2.16) with (2.2) that the functions 
J(s) and G{s) can be expressed in terms of the complex shear modulus /i(s) 
as 

As a consequence we are led to introduce an effective modulus defined as 

(2.19) 



1 



d 
di 



Git) + Gg 



Recalling that for perfect viscous fluid [Newton model) we have 

a{t) = 2r]oj^eit), (2.20) 

where tjq denotes the viscosity coefficient, similarly we are led to define, 
following Miiller [21] an effective viscosity as 

(2.21) 

where the dot denotes the derivative w.r.t. time t. 
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We easily recognize that for the Hooke model (2.4) we recover fi{t) = fiQ = 
m/2, and for the Newton model (2.5) ri{t) = rjo = bi/2. Furthermore 
we recognize that for a spring (Hooke model) the corresponding complex 
modulus is a constant, that is 

/2^(s)=/xo, (2.22) 
whereas for a dashpot (Newton model) 

/iAr(s) = r]oS = /ioSTo , (2.23) 

where tq = ?7o//io denotes a characteristic time related to viscosity. 
In order to avoid possible misunderstanding, we explicitly note hat the 
complex modulus is not the Laplace transform of the effective modulus but 
its Laplace transform multiplied by s. 

The appropriate form of the complex modulus /i(s) can be obtained recalling 
the combination rule for which creep compliances add in series while 
relaxation moduli add in parallel, as stated at the beginning of subsection 
2.2. Accordingly, for a serial combination of two viscoelastic models with 
individual complex moduli /ii(s) and /i2(s), we have 

^ ^ (2.24) 



/i(s) /ii(s) /22(s) 
whereas for a combination in parallel 

Jl{s) = -fli{s)+Ms)- (2.25) 



We close this section with a discussion about the definition of solid-like 
and fluid-like behaviour for viscoelastic materials. The matter is of course 
subjected to personal opinions. 

Generally one may define a fluid if it can creep indefinitely under constant 
stress (Je = oo), namely when it relaxes to zero under constant deformation 
{Ge = 0). According to this view, viscoelastic models of type II and IV are 
fluid-like whereas models of type I and III are solid-like. However, in his 
interesting book [33], Pipkin calls a solid if the integral of G{t) from zero to 
infinity diverges. This includes those cases in which the equilibrium modulus 
Ge is not zero. It also includes cases in which Gg = but the approach to 
the limit is not fast enough for integrability, for example G{t) ^ with 
0<a<last— )■ oo. 
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3 Fractional viscoelasticity 



The straightforward way to introduce fractional derivatives in hnear 
viscoelasticity is to replace in the constitutive equation (2.5) of the Newton 
model the first derivative with a fractional derivative of order v G (0,1), 
that, being e(O^) = 0, may be intended both in the Riemann-Liouville or 
Caputo sense. We refer the reader to Appendix A and Appendix B for the 
essential notions of fractional calculus for causal systems and for a discussion 
on initial conditions. 

Some people call the fractional model of the Newtonian dashpot (fractional 
dashpot) with the suggestive name pot. we prefer to refer such model to as 
Scott-Blair model, to honour to the scientist who already in the middle of 
the past century proposed such a constitutive equation to explain a material 
property that is intermediate between the elastic modulus (Hooke solid) and 
the coefficient of viscosity (Newton fluid), see e.g. [25] • Scott-Blair was surely 
a pioneer of the fractional calculus even if he did not provide a mathematical 
theory accepted by mathematicians of his time, as pointed out in [28] . 

It is known that the creep and relaxation power laws of the Scott-Blair 
model can be interpreted in terms of a continuous spectrum of retardation 
and relaxation times, respectively, see e.g. [27], p. 58. Starting from these 
continuous spectra, in [32] Papoulia et al. have interpreted the fractional 
dashpot by an infinite combination of Kelvin- Voigt or Maxwell elements 
in series or in parallel, respectively. We also note that Liu and Xu [21] 
have computed the relaxation and creep functions for higher-order fractional 
rheological models involving more parameters than in our analysis. 

3.1 Fractional derivatives in mechanical models 

The use of fractional calculus in linear viscoelasticity leads to generalizations 
of the classical mechanical models: the basic Newton element is substituted 
by the more general Scott-Blair element (of order u). In fact, we can 
construct the class of these generalized models from Hooke and Scott-Blair 
elements, disposed singly and in branches of two (in series or in parallel). 
Then, extending the procedures of the classical mechanical models (based on 
springs and dashpots), we will get the fractional operator equation (that is 
an operator equation with fractional derivatives) in the form which properly 
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generalizes (2.14), i.e. 



k=l 



a{t) 



m 



k=l 



dfk 



t{t), Vk = k + v-l. (3.1) 



so, as a generalization of (2.11), 



Jit) = J, + 5^ Jn {1 - h(t/r,,„)1} + J+fTl 



(1 + ^) 



Git) = G, + J2GnK [-(t/r.,„)1 + — 



(3.2) 



(1-z.) 



where all the coefficient are non-negative. Here F denotes the well known 
Gamma function, and E^, denotes the Mittag-Leffier function of order u 
discussed hereafter along with its generalized form in two parameters -E,^,^. Of 
course, for the fractional operator equation (3.1) the four cases summarized 
in Table 2.2 are expected to occur in analogy with the operator equation 
(2. 14). The definitions in the complex plane of the Mittag-Leffier functions 
in one and two parameters are provided by their Taylor powers series around 
z = 0, that is 



n=0 



n=0 



r(z/n + /i) 



(3.3) 



related between them by the following expressions, see e.g. [27^, Appendix 
E, 



E,{z) = E,Az) = 1 + zE,^^+,{z) , ^ E.^z") 



z--'E,Jz-). 



(3.4) 



In our case the Mittag-Leffier functions appearing in (3.2) are all of order 
V G (0.1] and argument real and negative, namely, they are of the type 

EA-it/m = ^(-1)" ^ ' y 0<u<l, r>0, (3.5) 

n=0 ' 

that for 1/ = 1 reduce to exp(— t/r). 
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Let us now outline some noteworthy properties of the Mittag-Leffler function 
(3.5) assuming for brevity r = 1. Since the asymptotic behaviours for small 
and large times are as following 

M-n ~ <; t y (3-6) 

t — )■ +00 - 



we recognize that for t > the Mittag-Lefiler function decays for short times 
like a stretched exponential and for large times with a negative power law. 
Furthermore, it turns out to be completely monotonia in < t < oo (that is 
its derivatives of successive order exhibit alternating signs like e~*), so it can 
be expressed in term of a continuous distribution of elementary relaxation 
processes: 

E.(-n = K.ir) dr . K.ir) = ± > ■ (3.7) 

For completeness we also recall its relations with the corresponding Mittag- 
Leffler function of two parameters as obtained from (3.4): 

E^i-t") = E,,,i-tn = 1 - r E,,,+,i-f) , ^ E,{-tn = -r-i E,,,(-r) . 

(3.8) 



We recognize that in fractional viscoelasticity governed by the operator 
equation (3.1) the corresponding material functions are obtained by using 
the combination rule valid for the classical mechanical models. Their deter- 
mination is made easy if we take into account the following correspondence 
principle between the classical and fractional mechanical models, outlined in 
1971 by Caputo and Mainardi [S]. Taking < u < 1 and denoting by r > a 
characteristic time related to viscosity, such correspondence principle can be 
formally stated by the following three equations where the Laplace transform 
pairs are outlined as well: 

itlrV 1 

5{t/r) - r ^ ^ _ ^srY , (3.9) 

1 1 {t/TY 11 , , 

t T ^ - — — ^ J"/ — — , 3.10 

' s sr r 1 + u) s {stY ' ^ ' 
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e-V^ ^ ^ EJ-(t/TY] - - -^^^L . (3.11) 

In the following, we will provide the creep compliance J{t), the relaxation 
modulus G{t) and the effective viscosity rj{t) for a set of fractional models 
which properly generalize with fractional derivatives the basic mechanical 
models discussed in Section 2, that is Kelvin- Voigt, Maxwell, Zcncr, anti- 
Zener, and Burgers, by using their connections with Hooke and Scott-Blair 
elements. Henceforth, for brevity we refer to all the basic models with their 
first letters, that is H, SB, KV , M, Z, AZ and B. Our analysis will be 
carried out by using the Laplace transform and the complex shear modulus 
of the elementary Hooke and Scott-Blair models. For this purpose we recall 
the constitutive equations for the H and SB models 

Hooke model : a{t)=me(t), Scott — Blair model 

and, setting 

^ 2 ' ^ 2// ' 
we get the complex shear moduli for these elements 

/ii^(6-) = /i , ^sb{s) = iAstY , (3.14) 

where u G (0, 1] and r > is a characteristic time of the SB element. So, 
whereas the H element is characterized by a unique parameter, its elastic 
modulus II, the SB element is characterized by a triplet of parameters, that 
is {ii,r,u}. 

3.2 Fractional Kelvin— Voigt model 

The constitutive equation for the fractional Kelvin- Voigt model (referred to 
as KV body) is obtained from (2.6a) in the form 

d^e 

fractional Kelvin — Voigt model : a{t) — me{t) + hi . (3.15) 

The mechanical analogue of the KV body is represented by a Hooke 
{H) element in parallel with a Scott-Blair {SB) element. The parallel 



-W = 6.-.(t), 

(3.12) 
(3.13) 
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Figure 4: Normalized creep compliance (a), relaxation modulus (b) and 
effective viscosity (c) for the KV body, for some values of the fractional 
power u in the range (0 < < 1), as a function of normalized time ^. The 
thick lines {u = 1) represent material functions and effective viscosity for the 
traditional Kelvin- Voingt body. 



combination rule (2.25) provides the complex modulus 

in 



2// ' 



(3.16) 
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where our time constant r reduces for u = 1 to the retardation time of the 
classical KV body, see (2.6b). Hence, substitution of (3.16) in Eqs. (2.18) 
gives 



Jkv{s) 



2fis 



ST 



1 + (st) 



Gkv{s) 



2/i 



[1 + [srY] . (3.17) 



By inverting the above Laplace transforms according to Eqs. (3.11) and 
(3.9), we get for t > 0, 



JMt) = ^ [1 - E.i-it/rY)] , GKvit) = 2/i 



1 + 



(3.18) 



For plotting purposes, it is convenient to introduce a non-dimensional time 
^ := t/r and define normalized, non-dimensional material functions. With 
J'^viO = f^JKvit)\t=T^ and G^y(^) = {1/ fi)G Kvit)\t=T£, these can be written 
in the non-dimensional form 



1 + 



(3.19) 

Finally, by a straightforward application of Eq. (2.21), recalling the 
derivative rule in (3.8) for the Mittag-Leffler function, the effective viscosity 
of the KV body turns out to be 



l-u 



(3.20) 



Figure H] shows plots of of JxyiO^ G'f^y{^) and rj'j^viO as a function 
of non-dimensional time ^. For z/ — )■ 1, the response of the fractional 
KV body reduces to that of a classical KV body. Taking into account 
that ^i(-0 
to J'KviO = 



the creep compliance reduces, in this limiting case, 
1 — e~^)/2. Recalling the Dirac's delta representation 
6(t) = t^^/T{0), the relaxation modulus reduces, for z/ — 1, to G'j^y{t) = 
2[1 + T6{t)] = 2[1 + 5(t/r)]. Finally, for the effective viscosity we obtain the 
classical exponential law Tj'j^viO = W^e note the effective viscosity exceeds 
the classical value since the early stage of creep. 
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3.3 Fractional Maxwell model 



The constitutive equation for the fractional Maxwell model (referred to as 
M body) is obtained from (2.7a) in the form 



d^a d^e 

fractional Maxwell model : cr{t) -\- Oi -— = hi -— . 

(JjV (JjV 



(3.21) 



The mechanical analogue of the M body is composed by a Hooke (H) 
element connected in series with a Scott-Blair (SB) element. Prom the 
series combination rule (2.24) we obtain the complex modulus as 



{srY 



1 



[ST 



l^ 



2ai 



2ii 



(3.22) 



where now our time constant r reduces for z/ = 1 to the relaxation time 
of the classical M body, see (2.7b). Hence, substitution of (3.22) into Eqs. 
(2.18) gives the Laplace transforms of the material functions 



Jm{s) 



1 



2 lis 



1 + 



Gm{s) 



2n {stY 
s 1 + {stY 



(3.23) 



By inverting the Laplace transforms according to Eqs. (3.10) and (3.11), we 
get for t > 0, 



1 + 



r(i + ^). 



GMit) = 2fx (- (t/r) 



(3.24) 



For plotting purposes, it is convenient to write non-dimensional forms of 

Juit) and Guit)- These can be obtained, following the example of the 
fractional KV body, by introducing a non-dimensional time with ^ = 
t/r and defining normalized, non-dimensional material functions J^iC) = 
lJiJM{t)\t=Ti and Gm(0 = (V/^) GM{t)\t=Ti, which provides 



1 + 



G'uiO-'^E^i-n , e=-- (3.25) 



Following this normalization scheme, the effective viscosity (in non- 
dimensional form) can be readily obtained from (2.21) and (3.25) as 



r(i + i^) ^1-. 



(3.26) 
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Figure 5: Normalized creep compliance (a), relaxation modulus (b) and 
effective viscosity (c) for the M body, for some rational values of u in the 
range (0 < < 1), as a function of normalized time ^. The thick lines, 
corresponding to z/ = 1, show the classical Maxwell response. 



Plots of J'm{C)) G^m(0 ri'^{C) are shown in Figure [5l Thick lines show 
material functions and effective viscosity in the limit of z/ i— )■ 1, i. e., when 
the response of the fractional M body degenerates into that of a classical 
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M body. In particular, from Eq. (3.25) we obtain J'm{C) = (1 + 0/2 and 
^'hiiO — 2e^^ and the effective viscosity is constant (^7^(0 — 1)) which 
denotes the lack of transient effects. For < u < 1, the effective viscosity 
always increases with time and, except during the very early stages of creep 
(^ fa 1) it exceeds the value corresponding to the classical M body {u — 1). 



3.4 Fractional Zener model 

The constitutive equation for the fractional Zener model (referred to as Z 
body) is obtained from Eq. (2.8a) in the form 



fractional Zener model : a{t) + ai 



(Fa 



m e{t) + h] 



'dF 



(3.27) 



The mechanical analogue of the fractional Z body is represented by a Hooke 
(H) element in series with a Kelvin- Voigt (KV) element. Here we indicate 
with III the shear modulus of the H body while with {ii2,T2,iy} the triplet 
of parameters characterizing the KV body. 

The material functions for the fractional Z body in the time domain can 
be derived following the same procedure outlined above for the fractional 
KVand M. However, the algebraic complexity increases because of the 
increased number of independent rheological parameters involved. The use 
of the combination rule (2.24), which holds for connections in series, provides 
the complex modulus 



1 



1 



1 1 

Hi /X2 1 + (sTaJ 



bi_ 
2ai 



m 



— , (3.28) 



where now our time constant T2 reduces for z/ = 1 to the retardation time 
Te of the classical Z body, see (2.8b). Hence, substitution of (3.28) in Eqs. 
(2.18) provides the Laplace transform of the creep compliance and relaxation 
modulus. For the creep compliance we have 



Jzis) 



1 

2^ 



1 1 



1 



1^2 1 + {sT2y 



(3.29) 



that can be easily inverted to obtain 



- + -{l-E,{-{t/T2) 
Pi A*2 



t > 0. 



(3.30) 
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Figure 6: Normalized creep compliance (a), relaxation modulus (b) and 
effective viscosity (c) for the Z body, for some values of the fractional power 
V in the range (0 < < 1), as a function of normalized time ^. Here we 
adopt the ratio = ^2/1^1 = 1- The thick lines {u = 1) represent material 
functions and effective viscosity for the traditional Zener body. 

For the relaxation modulus the following expression can be obtained 
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where 

= ^ , (3.32) 
and we have introduced an additional characteristic time 

< = YT7-'^^ '^ = ^- (3.33) 

Applying a partial fraction expansion to Eq. (3.31) and using Eqs. (3.8), the 
relaxation modulus in the time domain can be cast in the form 

Gz{t) = 2/i* (^^^ [E, (- it/TaY) + (t/r2)^E,,,+i (- {t/rS)] , (3.34) 

hence, after some rearrangement, we finally obtain 

Gz{t) = 2/i* [1 + r^E^i-it/TaY)] . (3.35) 

We now recognize that for u = 1 the constant reduces to the relaxation 
time To- (0 < To- < < oo) for the classical Z model, see (2.8b). 

Writing J^(0 = ft*Jz{t)\t=r2i and G'ziO = i^/ fJ'*)Gz{t)\t=r2^ the material 
functions for the Zener model can be written in the non-dimensional form 



G',{o = 2[^ + hiE,{-hn], ^ = -, 

(3.36) 

with hi = r^, /i2 = 1 + r^. 

As a consequence the effective viscosity of the Z body turns out to be 

Figure |6] shows plots of J'ziOy G'^iCj and ?7^(0 as a function of ^, for 
various values of v in the range < z/ < 1. In the limiting case of z/ = 1 
(thick curves), which corresponds to the classical Zener model, the material 
functions and the effective viscosity can be written in terms of exponential 
functions, in agreement with Eqs. (2.8b) above. In particular, since = 1, 
we obtain J'z{Cl = ^2(1 - V2 6"^) and G'ziO = 2(1 + 6"^^). Furthermore, 
the effective viscosity grows exponentially, since ?7^(0 = ^2 e^. We note 
that, for sufficiently long times (^ > 0.2 in Figure |6]), the value of ?7^(0 
for a fractional Zener model {0 < u < 1) always exceeds the classical Zener 
viscosity. 
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3.5 Fractional anti— Zener model 



The constitutive equation for fractional anti-Zener model (referred to as AZ 
body) is obtained from (2.9a) in the form 



fractional anti— Zener model : 



dt"" 



a{t) 



' d" ^ d^+^ 
^'df^^'dt^ 



e{t) 



(3.38) 

The fractional anti-Zener model results from the combination in series of 
a SB element (with material parameters {/x,ri,i/}) and a fractional KV 
element ({/x, r2, t^}). 

Using the series combination rule (2.24) we obtain 



1 



1 1 



1 + {ST2)\ 





f ai 


h2\ 












hi) 



1 



1 - r . , - = 2 I — - I , rr = T2 

(3.39) 

where now our time constant T2 reduces for v = 1 io the retardation time 
of the classical AZ body, see (2.9b). Using (2.18) the above equation easily 
provides the Laplace transform of the creep compliance 



Jaz{s) 



whose inverse is 



JAz{t) 



2/1 



2/xs 



+ 1 



1 + {ST2)\ 



;i - E^{-{t/T2r)) 



t > 0. 



(3.40) 



(3.41) 



The derivation of the relaxation modulus is more complicated but straight- 
forward. Starting from the general relationship (2.18) and using (3.40), by a 
partial fraction expansion we obtain 



Gaz{s) = ^ 
s 



(sroT + ( J 



2v 



[ST] 



where we have defined the new time constants 



1 + isfY 



riT2 



(3.42) 



(3.43) 



Hence, using (3.9) and (3.11), the inverse Laplace transform of (3.43) is 



GAz{t) = 2/i 



t>0. 



(3.44) 



hi 
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Figure 7: Normalized creep compliance (a), relaxation modulus (b) and 
effective viscosity (c) for the AZ body, for some values of the fractional power 
V in the range (0 < z/ < 1), as a function of normalized time ^. The particular 
case Ti = T2 is considered. The thick lines, pertaining to = 1, show results 
for the traditional AZ body. 

We now recognize that for v = \ the constant f reduces to the relaxation 
time Tfj (0 < < < oo) for the classical AZ model, see (2.9b). 
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From (3.41) and (3.44) it is easy to derive non dimensional expressions 
for the material functions in the time domain. The choice of the scaling for 
the time variable is of course arbitrary. Here we define 

t * riT2 

r* 



(3.45) 



n + T2 

With this choice, defining t/^^(0 = f^'JAz{t)\t=T'-£, and manipulating (3.41) 



we obtain 



where 



ic20' 



Cl 



C2 



(1 - (- (ciO^)) 

1 n 



(3.46) 



(3.47) 



1 + 1 + Tr T2 

Similarly, from (3.44), the relaxation modulus has the non dimensional form 



G'AziO = 2 



+ erE,i-iesO' 



(3.48) 



_r(l - u 

where G'^^(^) = (l//i) GAz{t)\t=T*£, and we have introduced the constants 



ei 



62 



(3.49) 



^ + rr ' ' {l + r^)t' " (l + r^)(l + r^)^ 

Using the normalization scheme above, the expression for the effective 
viscosity of the AZ body can be obtained from Eqs. (3.45) and (3.46) as 

1 



VAziO 



C2 (C2O 



Cl icO"'' E,,, (- icon 



(3.50) 



Figure [7] shows functions J'az{0^ ^!4z(0 and rj'^ziO as a function of ^, in the 
particular case Ti = T2 (r^ = 1 according to Eq. (3.46)). The behaviour for 
z/ = 1 (thick curves), which corresponds to the classical anti-Zener model, can 
be expressed in terms of elementary functions. Since r,- = 1, ci = C2 = 1/2, 
see Eq. (3.47). Hence, from (3.46), in this particular case the normalized 
creep function is J'az^O = 1/2 (l + + e~'^) , which, in the range of ^ 
values considered in Figure [TJ appears to be dominated by the linear term. 
In the same limit, since from Eqs. (3.49) we have ei = 1, 62 = 1/2, and 
63 = 1/4, omitting an additive 6(0 term, the normalized relaxation modulus 
decays exponentially with ^ according to G'^^iO = e~^^^/2. The effective 
viscosity, according to (3.50), turns out to be ?7^z(0 = 2/(1 + e~^/^), which 
approximately is increasing linearly in the range of C, values considered. 
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3.6 Fractional Burgers model 



The constitutive equation for the fractional Burgers model (referred to as S 
body) is obtained from (2.10a) as 



fractional Burgers model : 



1 + OiTr + 0,2-r-r— 



a{t)- 



■ d" , 

(3.51) 

The mechanical analogue of the fractional Burgers {B) model can be 
represented as the combination in series of a fractional KV element (with 
material parameters {/Xi,Ti,i/}) and a fractional M fractional element 

({/^2, 7-2,1^}). 

Using the expressions (3.16) and (3.22) for complex moduh of fractional KV 
and M bodies into the series combination rule (2.24) provides the complex 
modulus 



V'b{s) 



1 



1 + 



1 

H 

A«2 



1 



1 + {sT2)r 



(3.52) 



where the four constants //j and Tj are related in some way to the four 
coefficients Oj, hi of the constitutive equation with i — 1,2. 

Using (2.18), from (3.52) it is possible to obtain immediately the Laplace 
transformed creep compliance 



Jb{s) 



1 



1 1 
— + — 



+ 



1 1 



which can be easily inverted with the aid of (3.10) and (3.11), obtaining 



(3.53) 



JB{t) = - 



Ail A*2/ /xir(i + i/) 



-E^{-{t/T2Y) 



i > 0. 
(3.54) 



The computation of the Laplace-transformed relaxation modulus for the 
fractional B body requires more cumbersome algebra compared to Jb- Using 
(3.52) into (2.18) provides, after some algebra: 



Gb{s) = 2/xVo"s 



(3.55) 



Above, for the sake of convenience, we have introduced the variable 

z = {sToY , (3.56) 
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Figure 8: Normalized creep compliance (a), relaxation modulus (b) and 
effective viscosity (c) for the B body, for values of the fractional power 
1^ = 1, 0.75, 0.50 and 0.25 as a function of normalized time ^. The thick 
lines {ly — 1), show results for the traditional B body. 



where tq is defined as in the fractional AZ model (3.43), and several constants 
as follows: 

^* = ^ , (3.57) 
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1 + 1 + 1 + 1 + r^, 1 + < ' 1 + r^, (1 + r^^)^ " 



where = ri/T2 and = /Xi///2- 

Denoting by 2;i and ^2 the roots of algebraic equation 



(3.58) 



z^52 + + 5o = , (3.59) 

we obtain 

Gs(5) = 2/.X5--^ J]^, (3.60) 



where the evaluation of constants Gi — Gi{r^,rT) is straightforward (the 
details are left to the reader). Since from (3.58) it can be easily recognized 
that the discriminant of Eq. (3.59) is positive, the roots are real. Hence, 
observing that 6i > {i = 0,1,2), Descartes' rule of signs ensm^cs that 
< (i = 1,2). Setting Zi = —1/Pi with > and recalling (3.56), 
we obtain a convenient expression of the relaxation modulus in the Laplace 
domain: 

e.w-vE^.,|g^. (3.61) 

which can be easily inverted with the aid of (3.11), obtaining the relaxation 
modulus in the time domain as a combination of two Mittag-Leffler functions 
with different arguments: 

GB{t) = 2i^* J2 ^i^- (- (:^) ) ' (3-62) 

Useful non-dimensional forms for the material functions in the time domain 
can be easily obtained. For the creep compliance (3.54) we have 



4(0 = \ 



i + 4p|^-c;i?.(-(c,oi 



i ' * '''' 



r* Ti + T2 

(3.63) 

where Js(^) = JB{t)\t=T*(,, and 

4 = T^. 4=^. .. = ^. (3.64) 
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For the non-dimensional form for the relaxation modulus we obtain: 

G'niO = 2^^G.E^ {- (^J-^jy (3.65) 

where G'^{0 = iUP*) GBm=r'(, ^ = t/r*. 

By the above normalization scheme, it is possible to find the expression for 
the effective viscosity of the B body by using Eqs. (2.21) and (3.63) as 

v'siO = V ^ • (3-66) 

r(l + vf '^^'^ ^^'^^> + ^^^1 ^~ ^^^^^ ^ 

As an illustration of the behaviour of the material functions for the B body, 
in Figure E] we consider, for ri = 2x2 and /ii = fi2, J'siO^ ^b(0 and //^(O 
as a function of variable C,- While the results for the creep compliance (a) 
and the effective viscosity (c) are qualitatively similar to those obtained for 
the AZ model in Figure [TJ those pertaining to the relaxation modulus (c) 
differ significantly, being now removed the singularity for ^ i— )■ displayed in 
Figure [Tt for u ^ 1. Curves obtained for u = 1 (thick), which corresponds 
to the classical B model, can be expressed in terms of elementary functions. 
Since r^- = 2 and = 1, in this case we have Ci = 2/3, C2 = 1/3, c[ = = 
1/2. Hence, with these numerical values for the rheological parameters the 
creep function varies with time as 

J'AziO = ^ (l + I - ^ ^'^^^^) ■ (3-67) 

In the same limit, since Ei{—x) = e~^, the relaxation modulus G'^{^) reduces 
to the sum of two exponentially decaying functions, and its initial value is 
only determined by the numerical value of the elastic parameters fii and 
fi2- The effective viscosity, according to (3.66), for u = 1 asymptotically 
approaches a constant value for ^ — )■ oo as in the case of the AZ body; for 
all the other values of the fractional order rj'^iC,) increases indefinitely, as in 
the case of the M body. 



4 Conclusions 

Starting from the mechanical analog of the basic fractional models of 
linear viscoelasticity we have verified the consistency of the correspondence 
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principle. This principle allows one to formally get the material functions 
and the effective viscosity of the fractional models starting form the known 
expressions of the corresponding classical models by using the corresponding 
rules (3.9), (3.10), (3.11). 

In particular, in the present paper, we were able to plot all these functions 
versus a suitable time scale in order to visually show the effect of the 
order u G (0, 1] entering in our basic fractional models. This can help the 
researchers to guess which fractional model can better fit the experimental 
results. Presentation of the results has greatly benefited from a recently 
published Fortran code for computing the Mittag-Leffler function of complex 
argument [391 SO], coupled with an open source program for manipulating 
and visualizing data sets [H]. 

In Earth rheology, the concept of effective viscosity is often introduced to 
describe the behaviour of composite materials that exhibit both linear and a 
non-linear (i.e., non-Newtonian) stress-strain components of deformation), 
see e.g., Giunchi and Spada [13]). However, according to Miiller [31] and to 
the results outlined in this paper, it is clear that this concept has a role also 
within the context of linear viscoelasticity, both for describing the transient 
creep of classical models and for characterizing the mechanical behavior of 
fractional models. 

Transient rheological effects are currently the subject of investigation in the 
field of post -seismic deformations [U |29l [30], global isostatic deformations 
[36l [37] and regional sea level variations [38]. The fractional mechanical 
models described in this paper have the potential of better characterizing 
the time-dependence of these processes, also in view of the increased quality 
of the available geophysical and geodetic observations. 
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Appendix A. The two fractional derivatives in 
R+ 

For a sufficiently well-behaved function f{t) {t G R^) we may define the 
fractional derivative of order /3(m — l</3<m, mGN), see e.g. the 
lectures notes by Gorenflo and Mainardi [12] and the text by Podlubny [31] 
in two different senses, that we refer here as to Riemann-Liouville derivative 
and Caputo derivative, respectively. Both derivatives are related to the so- 
called Riemann-Liouville fractional integral of order a > 0, defined as 



rf{t) 



T(a) 



(Al) 



We note the convention = Id (Identity operator) and the semi-group 
property 



Ta+13 



a > 0, /3 > 0. 



{A.2) 



The fractional derivative of order /3 > in the Riemann-Liouville sense is 
defined as the operator which is the left inverse of the Riemann-Liouville 
integral of order /3 (in analogy with the ordinary derivative), that is 



D^I^ = Id, /3 > . 



(A.3) 



If m denotes the positive integer such that m — l</3<m,we recognize 
from Eqs. (A.2) and (A.3) f{t) := D"^ /'""^ f{t) , hence 



d^ 



df^ 



r(m-/3) Jo (t-r)/3+i-™ 
fit), 



dr 



m 



/3 



1 < /3 < m. 



m. 



For completion we define = Id. 

On the other hand, the fractional derivative of order /3 > in the Caputo 
sense is defined as the operator such that f{t) := J"*"^ /(t) , hence 



1 



D'. fit) 



m)/ 



r(m - P) Jq {t - r)/^+i-™ 
d"" 



dr , 



dp 



;fit) 



m 



1 < /3 <m, 



(A5) 



m 
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where f^"^^ denotes the derivative of order m of the function /. We point out 
the major utihty of the Caputo fractional derivative in treating initial-value 
problems for physical and engineering applications where initial conditions 
are usually expressed in terms of integer-order derivatives. This can be easily 
seen using the Laplace transformation, according to which 

m—1 

L{D^^f{ty,s} = s('f{s)-Y,s('-'-'f^'\0+), m-l</3<m, (A6) 

k=0 

where /(s) = L{/(t);s} = / q'^^ f{t)dt, s G C and := 
lim f^^\t). The corresponding rule for the Riemann-Liouville derivative 

i-i-0+ 

is more cumbersome: for m — l</3<mit reads 

m—1 

L {D^ fit)- s] = / m - [D' /^'"-''T /(0+) s"^-'-' , {A.7) 

k=0 

where, in analogy with (A. 6), the limit for t — O"*" is understood to be taken 
after the operations of fractional integration and derivation. As soon as all 
the limiting values f^^\0^) are finite and m — 1 < (3 < m, the formula (A.7) 
simplifies into 

L{DPfit);s} = s^fis). (A8) 

In the special case /W(0+) = for A; = 0,1 , m — 1, we recover the identity 
between the two fractional derivatives. The Laplace transform rule (A. 6) was 
practically the starting point of Caputo, see O |6], in defining his generalized 
derivative in the late sixties of the last century. 

For further reading on the theory and applications of fractional calculus we 
recommend the more recent treatise by Kilbas, Srivastava and Trujillo [20] . 
For the basic results on fractional integrability and differentiability we refer 
the interested reader to the survey paper by Li and Zhao |23] . 

Appendix B: Remark on initial conditions in 
classical and fractional operator equations 

We note that the initial conditions at t = 0+, a^^\0^) with h = 0,1, . . .p — 1 
and e^''\0~^) with k = 0,1, . . . q — 1, do not appear in the operator equation 
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but they are required to be compatible with the integral equations (2.1). In 
fact, since Eqs. (2.1) do not contain the initial conditions, some compatibility 
conditions at t = 0^ must be implicitly required both for stress and strain. 
In other words, the equivalence between the integral Eqs. (2.1) and the 
differential operator Eq. (2.14) implies that when we apply the Laplace 
transform to both sides of Eq. (2.14) the contributions from the initial 
conditions are vanishing or cancel in pair-balance. This can be easily checked 
for the simplest classical models described by Eqs. (2.6)-(2.9). It turns out 
that the Laplace transform of the corresponding constitutive equations does 
not contain any initial conditions: they are all hidden being zero or balanced 
between the RHS and LHS of the transformed equation. As simple examples 
let us consider the Kelvin- Voigt model for which p = 0, q = 1 and m > 0, 
see Eq. (2.6), and the Maxwell model for which p — q = 1 and m — 0, see 
Eq. (2.7). 

For the Kelvin-Voigt model we get sa{s) = me{s) + b]se{s) — e(0+)], so, for 
any causal stress and strain histories, it would be 

sJ(s) = — ^ <^ e(0+) = 0. (B.l) 

We note that the condition e(O^) = is surely satisfied for any reasonable 
stress history since = 0, but is not valid for any reasonable strain history; 
in fact, if we consider the relaxation test for which e{t) = Q{t) we have 
e(0''") = 1. This fact may be understood recalling that for the Voigt model 
we have Jg = and Gg — oo (due to the delta contribution in the relaxation 
modulus). 

For the Maxwell model we get 5(s) + a [sa{s) — o'(O^)] = h[se{s) — e(0^)], 
so, for any causal stress and strain histories it would be 

sJis) = ^ + 1 ^ a(7(0+) = 6e(0+) . {B.2) 

OS 

We now note that the condition aa{Q'^) = be{0'^) is surely satisfied for any 
causal history both in stress and in strain. This fact may be understood 
recaUing that for the Maxwell model we have Jg > and Gg — 1/ Jg > 0. 

Then we can generalize the above considerations stating that the compat- 
ibility relations of the initial conditions are valid for all the four types of 
viscoelasticity, as far as the creep representation is considered. When the 
relaxation representation is considered, caution is required for the types III 
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and IV, for which, for correctness, we would use the generahzed theory of 
integral transforms suitable just for dealing with generalized functions. 

Similarly with the operator equations of integer order we note that the 
initial conditions at t = O"*" for the stress and strain do not explicitly enter 
into the fractional operator equation (3.1). This means that the approach 
with the Caputo derivative, which requires in the Laplace domain the same 
initial conditions as the classical models is quite correct. However, if we 
assume the same initial conditions, the approach with the Riemann-Liouville 
derivative provides the same results since, in view of the corresponding 
Laplace transform rule (A. 8), the initial conditions do not appear in the 
Laplace domain. The equivalence of the two approaches was formerly 
outlined by Mainardi |26| and more recently noted by Bagley [1] for the 
fractional Zener model. 

We refer the reader to the paper by Heymans and Podubny [16] for 
the physical interpretation of initial conditions for fractional differential 
equations with Riemann-Liouville derivatives, especially in viscoelasticity. 
In such field, however, we prefer to adopt the Caputo derivative since it 
requires the same initial conditions as in the classical cases. By the way, for 
a physical view point, these initial conditions are more accessible than those 
required in the more general Riemann-Liouville approach, see (A. 7). 
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